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I. INTRODUCTION 



Experiments on sliding drops on solid substrates or on moving contact lines show that the driv- 
ing force must exceed a nonzero threshold in order that motion results Jl-4]. These observations 
are at variance with the theoretical predictions obtained for ideally smooth and homogeneous 
substrates for which motion starts for arbitrarily small driving force [Q, SJ9]. This discrepancy 
between experiment and theory is believed to be caused by chemical heterogeneities and/or topo- 



graphic roughness that are always present on real substrates lUOj: 



ityO, 



4fl. A 



inite force is then needed 



12L113 



15 



m. 



to depin the contact line or entire drop from the heterogeneity 

More generally, micro- or mesoscale heterogeneities are expected to affect the macroscopic 
movement of drops. For instance, they are responsible for contact angle hysteresis [[J AH, the 
roug hening of contact lines II ll l2l 120-42211 . and the stick- slip motion of weakly driven contact lines 



[17, 



23|]. Heterogeneities 



patterns in dewetting 11241 - 12611 



may also slow down or stop entirely the long-time coarsening of drop 



These considerations apply to drops or liquid fronts (contact lines) moving down an inclined 
substrate due to a gradient of potential energy [3j,|5j, |27|,|28|]. However, they also apply to motion 
of drops or fronts caused by a temperature gradient along the substrate (resu lting in Marangoni 



forces due to surface tension gradients) 11291. 13011. or by wettability gradients [8, 



10, 



31 -[ 



3IH34J]. 



One approach to real (i.e., non-idealized) substrates is to consider the limit of random hetero- 



geneities II20M22| 



fects 111, 



35 



37, 



35 



360 . Another approach focuses on the effect of individual well-defined de- 



38]. Recently, the latter approach was extended to study the depinning dynamics 
of drops on substrates with a periodic array of precisely specified defects J39-41]. The approach 
employs a thin film evolution equation with a spatially modulated disjoining pressure and enables 
one to (i) study the depinning transition employing tools from dynamical systems theory and bifur- 
cation theory, and (ii) investigate the dynamics of the stick-slip motion that occurs after depinning 
on substrates with many defects. Contact line and drop motion on regularly patterned substrates 
i42l445|l is in fact of considerable importance in various micro-fluidic and nanotechnological ap- 
plications ||46l- l50|l . Single-cell assays, i.e., miniaturized devices for cell biology that consist of 
chemically and/or physically structured substrates, provide a good example. In this device each 
cell may be confined in an individual 'reaction chamber', e.g., a drop of water on a hydrophilic 
spot. A pattern of such spots allows for parallel analysis of a large number of cells 115 ill . Driven 
drop motion on a regular heterogeneous substrate with a well-defined wettability period is also 



2 



related to the motion of drops of partially wetting liquid on a horizontal rotating cylinder 113211 . 

At first sight the qualitative behavior of driven drops on substrates with a well-defined array 
of defects described in previous studies of the problem in a two-dimensional (2d) setting 1139 . 
I40I1 agrees well with the results of three-dimensional (3d) computations [41]. In a 2d system 
pinned drops can depin by two different mechanisms depending on the wettability properties of 
the defect, the drop size, and the driving force. The drops are either pinned by a hydrophilic defect 
at their back or by a hydrophobic defect at their front. For the parameters investigated in Ref. fl39h 
the following depinning scenarios are found. In the hydrophilic case the pinned drop stretches 
quasistatically as the driving force increases but eventually loses stability through a 'sniper' (or 
Saddle-Node Infinite PERiod) bifurcation resulting in depinning. For forces larger than this critical 
force the drop slides in a periodic motion over the periodic array of defects. Theory implies that the 
mean drop speed beyond depinning should vanish as the square root of the distance from the sniper 
bifurcation 115 311 and this behavior is indeed observed in simulations. Each period of the resulting 
stick-slip motion consists of two distinct evolution phases that take place on distinct timescales: 
the drop slowly stretches away from the defect but once it breaks away it slides rapidly over to 
the next defect. The time scales for the 'sticking' and 'sliding' phases differ greatly close to the 
bifurcation, and the behavior that results resembles closely the experimentally observed stick-slip 
motion. The situation is richer for hydrophobic defects that pin the drop by blocking it in front. 
In this case, in addition to the steady state sniper bifurcation, depinning can also occur via a Hopf 
bifurcation depending on the details of the defect Q39Q. 

A recent study of the depinning of 3d drops from hydrophobic and hydrophilic line defects j4l\ \ 
employing continuation and time-stepping algorithms [|54f | establishes a qualitative similarity be- 
tween the 2d and 3d cases and supports the widely held expectation that studies of 2d thin film 
systems provide useful information about more realistic 3d systems. However, significant differ- 
ences do exist. These are mostly related to the additional degrees of freedom of the 3d system. 
For instance, the 3d drop can, under appropriate circumstances, employ depinning pathways via 
morphological changes that are unavailable to a 2d drop. In the case of a hydrophilic defect the 
backward thread that connects the drop to the defect may gradually thin whereas in the hydropho- 
bic case the drop may 'probe' the barrier locally by sending out an advancing protrusion over the 
defect 114 ill . In fact, the 3d system has a number of features that have no counterpart in the 2d 
system, beyond the details of the depinning behavior of individual drops. The present study is 
dedicated to their analysis. 
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Our starting point is the observation that all 2d drop solutions correspond in a 3d setting to 
spanwise-invariant ridge solutions. It follows therefore that all results of I39I I40I1 remain valid 
in a 3d setting provided one imposes translation symmetry in the spanwise direction. Thus stable 
ridge and drop solutions of identical liquid volume may coexist. Either of these solution types may 
become unstable as parameters are varied giving rise to branches of solutions that have not yet been 
studied. In the limit of zero driving (e.g., a horizontal substrate when gravity is the driving force) 
and a homogeneous substrate this question is related to the characteristics of the Plateau-Rayleigh 



580. 



instability II55L 15611 of a static ridge (sometimes called a truncated cylinder) H57L 

In a driven 3d system on a heterogeneous substrate the question translates into (i) an investiga- 
tion of the lateral stability of pinned and sliding spanwise-invariant ridges, and (ii) the longitudinal 
stability of streamwise rivulets whose diameter is modulated by the heterogeneities. In both cases 
we expect changes of stability that give rise to 'new' branches of (i) spanwise modulated ridges 
and (ii) 'wavy rivulets', respectively. These states have no counterpart in the 2d case. We are 



interested in particular in the relation between the stability of the rivulets [59] and the existence of 



depinned sliding drop solutions. The combination of these aspects of the problem together with 2d 



and 3d depinning characteristics obtained previously |39J-|41|] provides a fairly complete picture 
of the 3d problem and of the relation among static and sliding ridges, static and sliding drops and 
steady or wavy rivulets. 

The paper is organized as follows. In section |n] we summarize the thin film model we use for 
the study of drop depinning in 3d and introduce the numerical tools we employ. In section [Till 
we discuss the results, first for drops and ridges on a horizontal substrate (section UlI A|) . and then 
for ridges on an inclined substrate (section IIII B|) . followed by drops and drop-like states on an 
inclined substrate (section 1111 CI) . We interpret our results in section [TV] followed by conclusions 
in section |V] 

H. MODEL AND NUMERICAL TOOLS 
A. Lubrication equation 

The partial differential equation governing the time evolution of the profile of a thin liquid film 



was discussed in depth in the 2d case in Ref. [39] and adapted to the 3d case in Ref. [41]. Briefly, 



we consider a layer or drop of partially wetting liquid (with a small equilibrium contact angle) on a 
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FIG. 1: (color online) Sketch of the geometry of the problem: A ridge or drop is pinned at a stripe-like 
hydrophobic defect, i.e., a heterogeneous wettability that depends only on the x-coordinate. A driving force 
\i acts in the x-direction. 

fiat chemically inhomogeneous two-dimensional solid substrate (see sketch in Fig.Q]). Long-wave 
theory allows us to derive an evolution equation for the film thickness profile h(x, y, t) directly 



from the Navier-Stokes and continuity equations [|60l. 16 We use no-slip and no penetration 
boundary conditions at the substrate, and the equilibrium of tangential and normal stresses at the 
free surface. The wettability properties are incorporated as a disjoining pressure D21, 160] . In the 
presence of a small (for consistency with the long-wave approximation) lateral body force in the 
^-direction we obtain the non-dimensional equation 



d t h = -V • {Q(h) [V (Ah + Il(h, x)) + yue x ]} 



(1) 



where V = (d x ,d y ) and A = d 2 xx + dy y are the planar gradient and Laplace operator, re- 
spectively, with (x, y) denoting the downstream and spanwise directions. The mobility function 
Q(h) = h 3 corresponds to a parabolic velocity profile (Poiseuille flow). Capillarity is represented 
by Ah (Laplace pressure). The position-dependent wetting properties are incorporated via a y- 
independent disjoining pressure H(h, x) in order to focus on stripe-like defects. Of the different 
functional forms for H(h) found in the literature |2L[62|], many allow for the presence of a precursor 
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film of thickness 1-10 nm on a 'dry' substrate and these are used to describe partial wetting. In this 
way 'true' film rupture in dewetting and the stress singularity at a moving contact line are avoided. 
Here, we employ long-range apolar van der Waals interactions combined with a short-range polar 
electrostatic or entropic interaction J2, 63I 64], leading to the dimensionless disjoining pressure 



n(/i,x) = A-[l + e£(x)]e-\ (2) 

where the parameter b > indicates the relative importance of the two antagonistic interactions. In 
the following we employ this particular form of the disjoining pressure with the parameter b fixed 



1. |3^-41] b 



that any qualitatively similar disjoining pressure II vie 



ds like results, as becomes clear when 



comparing, for example, the dewetting results in I164l-l69n. In Eq. ©, e and are the strength 
of the wettability contrast and the shape function of the heterogeneity, respectively. For the latter 
we use a Jacobi elliptic function to model periodic arrays of localized defects, specifically 

= 2{cn[2K(k)x/L x , k}} 2 - C, (3) 

where k is the modulus of the elliptic function and K{k) is the complete elliptic integral of 
the first kind. As k — > 1 Eq. © describes localized defects. Throughout the paper we take 
k = (1 - KT 6 ) 1 / 2 with C = 0.241136 so that £(x)dx = 0, i.e., so that the mean of the dis- 
joining pressure does not depend on the wettability contrast e. In addition, we take the period L x 
sufficiently large to avoid interactions between adjacent drops/defects. The resulting wettability 
profile is shown in Fig. Ob), lower panel. The wettability contrast e > [e < 0] indicates 
a hydrophobic [hydrophilic] defect, i.e., it tells us whether the defect is less [more] wettable than 
the surrounding substrate. Thus H(h, x) allows us to incorporate a stripe-like wettability pattern in 
the theory and to study the influence of chemical substrate heterogeneities or defects via a spatial 
modulation of the material parameters involved. This variation must take place on length scales 
much larger than the film thickness for consistency with the long-wave approximation QVOQ. 

In nondimensionalizing the model to arrive at Eq. (00) we have used the length scales I for 
the film thickness and 'k for the (2, y) coordinates, the time scale 3r)~//K 2 l for the time and 
the force scale ^/Ty/n 3 / 2 for the force. The length I corresponds to a characteristic scale for 
the thickness of the wetting layer, while 7 and rj are the surface tension and viscosity of the 
liquid, respectively; k is a typical energy density scale related to wettability. The ratio of the 



vertical and horizontal length scales used, ^Jln/'y, corresponds to the smallness parameter in the 



lubrication approximation. This parameter is also closely related to the equilibrium contact angle 
on a homogeneous substrate, cf. Q8Q. In writing Eq. £[]) we assumed that the lateral driving force, of 
dimensionless strength fi, does not depend on the film thickness. This is the case for gravitational 
or centrifugal forces. For driving forces that depend on the film thickness the mobility for the force 
will differ from that for the pressure term. This is the case for a purely lateral temperature gradient 
or a gradient of an electric field for a dielectric liquid in a capacitor. Note that the long-wave 
scaling used here implies that the dimensionless contact angle and driving force /i may be of order 
one. 

In the literature one finds two different ways of counting spatial dimensions in the problem 
at hand. On the one hand, focusing on the mathematical structure of Eq. (QQ) one distinguishes 
between one-dimensional (h depends on x only) and two-dimensional (h depends on x and y) 
cases. On the other hand, one may count the physical dimensions and refer to the situation where 
the film thickness depends only on x [depends on (x, y)] as the two-dimensional (2d) case [three- 
dimensional (3d) case]. Here we follow the latter convention. 



B. Numerical schemes and parameters 



Based on Eq. (OQ) and the disjoining pressure © the 3d depinning behavior is analysed as 
follows: Steady-state solutions (pinned drops) and their stability are determined using continua- 
tion techniques and the stick-slip motion beyond the depinning threshold is analysed using time- 
stepping algorithms. Note that in the 2d case B9L M-Qfl an explicit scheme suffices for the latter and 
continuation can be performed using the package AUTO flV ll — iV3h . This is not possible in the 3d 
case where an accurate and effective time simulation of Eq. (OQ) remains a challenge, leading to a 
number of numerical issues lr74H76ll . For these reasons we employ here an exponential propaga- 



tion approach based on the exact solution of the linearized equation at each time- step QVVU - The 
required exponentiation of the Jacobian matrix is performed using a Krylov reduction based on 
the Cayley-Arnoldi algorithm II54I1 . The approach allows for a very good estimate of the optimal 
timestep in both the slow and fast dynamical regimes. This is of paramount importance since 
close to the depinning transition typical time scales vary over many orders of magnitude. The 
same Cayley-Arnoldi algorithm is employed in our tangent predictor-secant corrector scheme for 
the continuation of steady 3d drop states Q54Q . This approach is advantageous as it allows us to 



perform all tasks arising in a bifurcation analysis simultaneously. This includes the computation 
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of the kernel of the Jacobian to detect bifurcations and the stability analysis of the steady states. 
Steady solutions are characterized by their L 2 norm 



||M|| 



_L X Ly JQ 



(h(x, y) — h) 2 dxdy 



1/2 



(4) 



while time-dependent states such as sliding drops or ridges are characterized by their temporal 
period T and time-averaged norm 



||M|| 



T L x L y jo jQ 



(h(x, y, t) — h) 2 dxdydt 



1/2 



(5) 



In this paper we focus on a regular array of hydrophobic stripe-like defects with a wettability 
contrast e = 0.3 with periodic boundary conditions in the downstream direction with period L x = 
40 containing a single defect. Periodic boundary conditions are employed in the spanwise direction 
as well, with period L y that may be varied to study drops of different volume and spanwise ridges 
of different length. In these circumstances L y is used as a bifurcation parameter and is then denoted 
by the symbol L. In this paper the values of L are restricted by the requirement that the domain 
contains no more than one drop in the spanwise direction. The mean film height is fixed at h = 1.2, 
compared with the precursor film height h — 1. 

III. RESULTS 

A. Drops and ridges on a horizontal substrate: \i = 



On a heterogeneous substrate with a periodic array of hydrophobic defects the unique stable 
solution in the 2d case corresponds to a 2d drop sitting in the middle between two hydrophobic 
defects |39rj4jJ]. In the 3d case this solution corresponds to a ridge with translational invariance in 
the spanwise direction. Such a ridge may, however, be unstable to a Plateau-Rayleigh instability if 
the spanwise system size exceeds a critical value L c as occurs for ridges on homogenous substrates 



B57L 



58|] . For cylindrical liquid bridges between two solids the primary bifurcation is an imperfect 



subcritical pitchfork whose details depend on the particular setting (with/without gravity, fhermo- 
capillarity etc.), the mode number, the boundary conditions at the two supports and their geometry 

£-82] 



26. 



57. 



[|78L l79|1 . Less is known for a ridge on a solid substrate, where most results [|25 . 
concern linear stability. For example, ridges on striped heterogeneous substrates can be stabilized 
w.r.t. the Plateau-Rayleigh instability for any L y if they sit on a hydrophilic stripe or between hy 



drophobic stripes of sufficient wettability contrast [26J]. In the nonlinear regime we expect to find a 



(subcritical) pitchfork of revolution (due to translation invariance in the spanwise direction in our 
periodic setting) when the spanwise system size L = L y increases. Bulge states with large contact 
angles are studied in I 83L 
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FIG. 2: (color online) (a) Steady ridge (SR) solutions on a horizontal substrate with hydrophobic line 
defects as a function of the wettability contrast e. Solid [dashed] lines indicate stable [unstable] solutions. 
The symbols ± indicate the stability of the branches w.r.t. 2d perturbations, with — indicating stability and 
the +'s indicating the number of unstable eigenmodes. (b) Profiles of the solutions SRi, SR2 and SR3 at 
e = 0.3 (top panel) together with the heterogeneity profile ^(x) (lower panel). Here and in subsequent 
figures the hydrophobic defect is always centered at x = (x = L). Parameters: h = 1.2, L x = 40. 



Figure [2t a) shows the 2d solutions on a horizontal substrate (/i = 0) as a function of the wet- 
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FIG. 3: (color online) Bifurcation diagrams for (a) e = and (b) e = 0.3 showing the loss of stability of the 
SRi solutions (horizontal red line) with respect to 3d perturbations when the spanwise period L increases 
together with the resulting branch of 3d states labeled SDi (black line). Solid [dashed] lines indicate stable 
[unstable] solutions. Solution profiles at locations indicated by open squares in panel (b) are shown in 
Fig.lH Parameters: fi = 0,h = 1.2, L x = 40. 

tability contrast e and indicates the presence of three 2d states when e = 0.3. The solutions take 
the form of steady ridges (SR) of different L 2 norm and differ in both the location of the ridge rel- 
ative to the hydrophobic heterogeneity and their linear stability properties. The largest amplitude 
solution, labeled SRi, consists of a ridge confined midway between adjacent hydrophobic defects 
[Fig. Etb)]; this solution is stable with respect to 2d perturbations. In contrast, the solutions SR2 
and SR 3 are unstable. Of these, SR 2 consists of ridges superposed on top of the defects, a con- 
figuration that is expected to be unstable. The lowest amplitude solution SR 3 is characterized by 
minima in the film profile h(x) at both the defects and half way between them [Fig. (2b)]. This 
solution is also unstable. 

Figure [3] extends these results to 3d structures on a horizontal substrate (fi = 0). To obtain 
the figure we start with a stable SRi solution in Fig. |2] with either e = (homogeneous substrate) 
or e = 0.3 (hydrophobic line defect) and study its linear stability with respect to 3d volume- 
conserving perturbations with period L in the spanwise direction. A symmetry-breaking Plateau- 
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FIG. 4: (color online) Steady drop-like states SDi at locations indicated by open squares in Fig. [3jb) in 
terms of contours of constant h(x,y). (a) The saddle-node bifurcation at L® nl = 24.44, (b) the unstable 
branch at L = 25.09, (c) the stable branch at L = 25.08 and (d) L = 30.26. The downslope direction x is 
from top to bottom, with y horizontal. Parameters: fi = 0, e = 0.3, h = 1.2, L x = 40. 

Rayleigh-like instability sets in for L > L° cl (the superscript zero refers to fi = 0, i.e., the absence 
of a driving force). At L° cl a branch of steady 3d solutions bifurcates from the 2d states when an 
eigenvalue of double multiplicity becomes unstable. The double multiplicity is a consequence of 
0(2) symmetry of the SR stability problem under translations in y modulo the period L together 
with the reflection y — > —y and is not indicated in Fig.[3]or subsequent figures. The 3d solutions 
that result take the form of unstable ridges modulated in the spanwise direction [see Fig.HJb)] and 
are unstable. With increasing modulation amplitude these solutions turn around at a saddle-node 
bifurcation (at L = L® nl ) and acquire stability. At approximately the same location the solution 
ceases to resemble a spatially modulated ridge and begins to resemble a steady drop-like (SD) 
state. In the following we use the notation SD to refer to all 3d states, even near the bifurcation 
L° cl , where their appearance is ridge-like. It follows that for L < L° nl the only stable solution is 
the 2d ridge state SRi that is invariant w.r.t. spanwise translations. For L > LP cl the ridge is linearly 
unstable and decays into the drops SDi that constitute the only stable solution in this regime. In 
between, i.e., in the range L® nl < L < L Q cl , both the SRi and the larger amplitude SDx solutions 
are linearly stable, while the unstable subcritical SDi branch of modulated ridges corresponds to 
unstable threshold solutions separating the two stable solutions. Comparison of Figs.[3ta,b) shows 
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FIG. 5: (color online) Bifurcation diagram showing secondary bifurcations to 3d states from the branches 
SRi, SR2 and SR3 (horizontal red lines) of 2d steady ridges, and the associated branches of 3d steady 
drop-like states SDi, SD2, SD3 and SD4 (heavy black lines). Solid [dashed] lines indicate stable [unstable] 
solutions. Solution profiles at locations indicated by open squares are shown in Fig. [6j Parameters: fi = 0, 
e = 0.3, h = 1.2, L x = 40. 

that presence of the defect shifts the bifurcation to 3d states and the saddle-node bifurcation on 
the resulting branch of 3d states to larger values of the parameter L than required for the Plateau- 
Rayleigh instability on a homogeneous substrate. In fact, this stabilizing effect is much more 
pronounced in the hydrophilic case (e < 0, not shown). 

Figure |4] shows a sequence of solutions along the SDi branch of 3d states created at L° cl when 
/i = 0, corresponding to locations indicated by open squares in Fig.^b), i.e., for e = 0.3. Fig- 
ure BJa) shows the solution at the saddle-node (L® nl = 24.44), with (b) showing the unstable 
ridge-like state below the saddle-node (L = 25 < L° cl = 26.6) and (c) showing the corresponding 
stable drop-like state above the saddle-node, at the same value of L. FigureHfd) shows the (stable) 



12 



1.5 



2.5 



3.5 




(a) 



(b) 



(c) 



(d) 



(e) 



FIG. 6: (color online) Steady drop-like states at locations indicated by open squares in Fig. [5] in terms of 
contours of constant h(x, y). (a) SD 2 at L = 29.90, (b) SD 3 at L = 30.11, (c) SD 4 at L = 19.64 (close to 
SR2), (d) SD4 at L = 17.17 (close to SR3), and (e) SD5 at L = 30.25. The downslope direction x is from 
top to bottom, with y horizontal. Parameters: /j, = 0, e = 0.3, h = 1.2, L x = 40. 



drop state at L = 30 > LP cl . We emphasize that the notion of stability is limited to linear stability 
with respect to perturbations with spatial period L. It turns out that the large amplitude drop branch 
is typically unstable w.r.t. coarsening, i.e., to perturbations with periods that are integer multiples 
of L. However, we do not pursue questions related to coarsening in the present work. 

Figure [5] shows, in addition to the SDx states, four additional branches of drop-like states. Of 
these the states labeled SD 2 and SD 3 lie on a pair of unstable branches that are connected via a 
common saddle-node bifurcation at L sn2 ps 25.0 but that are disconnected from the ridge states 
SRi, SR2 and SR3 (Fig. |5]). The states labeled SD 2 consist of drops sitting on the defect with a 
maximum on the defect while SD 3 consist of drops sitting on the defect with maxima on either 
side of the defect [Fig. [6£a,b)], just as for the corresponding SR 2 and SR 3 states. In addition, we 
find two further branches of drop-like states, hereafter SD 4 and SD 5 , both created as a result of 
symmetry-breaking bifurcations of the SR 2 and SR 3 states. Following 118411 we refer to these states 
as varicose and zigzag states. The former are shown in Figs.[6tc,d) at two locations along the SD 4 
branch; the latter solution is shown in Fig. |6te). As shown in Fig. |5] the varicose [zigzag] branch 
bifurcates subcritically from SR 2 at L° 2 = 19.95 [L° 4 = 34.52] and connects to SR 3 supercritically 
at L° 3 = 16.05 [L° c5 = 23.19]. Note that the SD 4 branch is twice unstable while the SD 5 branch 



is four times unstable. As explained in 118411 this is a consequence of the fact that the zigzag state 
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is unstable with respect to two different varicose modes. Technically, the unstable eigenvalue of 
double multiplicity splits into two eigenvalues when moving from the SR branch onto the SD 5 
branch. 



B. Ridge on an inclined substrate: fi > 

Once the driving force [i becomes nonzero, the ridge or drop solutions will be displaced from 
their symmetric location midway between a pair of hydrophobic defects. This is a consequence 
of the term —^id x Q(h) in Eq. (Q~|) that breaks the reflection symmetry in the x -direction. On a 
homogeneous substrate the solutions acquire an asymmetric shape w.r.t. x — > —x and slide down- 



stream like a solitary wave of constant shape |5il85|]. On a heterogeneous substrate, however, the 



wettability defect may pin such a sliding solution depending on the relative size of the parameters 
/i and e. In the present case of a hydrophobic defect with e = 0.3 the ridge or drop will be blocked 
at its front. A finite force fi > /idcpin is needed to overcome the pinning effect of the defect in order 
that the ridge/drop may move. Beyond this depinning transition ridges/drops slide along the plane 
but do so with a nonuniform speed since both their shape and speed are modulated periodically as 



they pass individual defects in the periodic defect array 11391 - 14 111 . In 2d the depinning itself is either 
related to a saddle-node infinite period (sniper) bifurcation or to a Hopf bifurcation. As already 
mentioned the sniper bifurcation generates stick-slip motion just above /idepin since the ridge sticks 



for a long time to a defect before sliding to the next one 114011 . 

In this section we seek to elucidate the depinning process for 3d drops and relate it to the 
corresponding process in 2d. We begin with stability and depinning of the 2d states SR. 



1. Depinning in two dimensions 

Figure [7ja) presents the bifurcation diagram for e = 0.3 (hydrophobic defects) as a function 
of the driving force /i obtained using numerical continuation and direct numerical simulation in 
time. The figure shows 2d states only, with the L 2 norm © for steady ridges (SR) and the time- 
averaged L 2 norm © for the stick-slip states (SSR). As soon as /i ^ the profiles of the SR states 
become markedly asymmetric as shown for SRi in Fig. He). The figure illustrates the dramatic 
shift in the position of the dominant maximum of the solution towards the downstream defect 
arising from the competition between the driving force and the blocking effect of the defect that 
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FIG. 7: (color online) (a) Bifurcation diagram showing the amplitude of the three steady 2d states identified 
in Fig. |2] as a function of fi. Solid [dashed] lines indicate solutions that are stable [unstable] with respect to 
perturbations with a real eigenvalue. At larger [i only SRi exits, (b) Zoom of (a) showing in addition the 
branch of stable stick-slip states (SSR) (open triangles); the inset shows the time taken to travel between 
successive defects, also as a function of /i. The open circles denotes solutions that are stable with respect to 
oscillations. The symbols ± indicate the stability of the branches w.r.t. 2d perturbations, with — indicating 
stability and the +'s indicating the number of unstable eigenmodes. (c) Solution profiles along the upper 
SRi branch showing the blocking of the ridge by the downstream defect (the driving force acts towards 
the right), (d) The real part of the leading eigenvalues of the SRi states as a function of /i. The solid 
[dashed] lines correspond to real [complex] eigenvalues. The three different line widths correspond to the 
three sub-branches separated by saddle-node bifurcations. Parameters: e = 0.3, h = 1.2, L x = 40. 

prevents downstream motion. As the driving force fi increases the ridge profile steepens (i.e., the 
norm © increases). However, close to the saddle-node bifurcation at [i = fj%ni ~ 0.01157 the 
norm starts to decrease, and for the parameters used the saddle-node bifurcation is preceded by a 
Hopf bifurcation at /U^opf ~ 0.01154 [Fig. He)]. This bifurcation leads to a subcritical branch of 
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unstable small amplitude time-periodic "rocking" states hereafter referred to as oscillating ridges 
(OR). We conjecture that with decreasing \i these states undergo a global bifurcation to translating 
states at some ji = fi 2d (fi 2d < /i^pf < /^sni) These states differ, however, from the usual stick- 
slip states generated via the sniper bifurcation since the ridge must spend considerable time in the 
rocking state before shifting rapidly to the next defect downstream. 




(c) x x (d) 

FIG. 8: Space-time plots of the time evolution of stick-slipping ridge (SSR) solutions beyond depinning via 
a Hopf bifurcation over one period in space and time at (a) \i = 0.011414 close to the saddle-node of the 
SSR branch at = 0.01141303 with temporal period T = 2000.6, (b) \i = 0.012 with T = 875.0, (c) 
H = 0.02 with T = 345.0, and (d) far from depinning, /i = 0.06 with T = 116.8. 

Figure |7£b) zooms into the region of the bifurcation diagram where stable stick- slip ridges 
(SSR) are present, while Figs. [8] and [9] show several of these states as a function of time over one 
period T in terms of space-time plots (Fig. [8]) and L 2 norm (Fig.©. Here the 'period' refers to the 
time required for the ridge to slide from one defect to the next. Figure Hb, inset) suggests that the 
period T diverges logarithmically as fi — > fi g as expected from a global bifurcation. Examination 
of the SSR profile near /i = fi g when it is almost pinned suggests that the SSR approaches close 
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FIG. 9: (color online) The L 2 norm as a function of time for selected solutions on the branch of stick-slip 
ridges (SSR). Time is scaled by the period T listed in the legend together with the corresponding values of 

to an SRi on the intermediate segment of the SRi branch (not shown). The leading eigenvalues 
along the SRi branch are shown in Fig. |7£d) with the heavy red lines indicating the eigenvalues 
along the intermediate segment. For fi fa fi g the states on this segment have a single unstable 
eigenvalue \ u ~ 0.0101 and a leading stable eigenvalue X s « —0.00054. Thus \ u + \ s > 
and under these conditions standard theory shows that any global bifurcation involving states on 
the intermediate segment of the SRi branch must involve unstable periodic states (see the 3d case 
below for further discussion). This conclusion is consistent with the behavior of the SSR norm 
[Fig- Ob)] which suggests that the SSR states do indeed undergo a saddle-node bifurcation prior 
to any global bifurcation and so are unstable near [i = \i g . 

The leading eigenvalues on the intermediate segment of SRi are created by the collision on 
the positive real axis of the unstable complex eigenvalues responsible for the Hopf bifurcation on 
the upper segment of the SRi branch. After this collision the eigenvalues remain real with one 
continuing to increase as one follows the branch while the other crosses back into the negative 
half plane at the saddle-node /i^. At this point the steady branch turns back towards smaller [i 
[Fig. 13b)]. The bifurcation at /i^i is a standard saddle-node bifurcation and not a sniper. This 
is because for the parameter values used the depinning transition involves the unstable OR states 



0.3 
0.2 



17 



instead of the SRi states, i.e., with the appearance of the Hopf bifurcation the depinning transition 
moves from the steady states SRi to the time-dependent OR states. At a further saddle-node at 
/i^ 2 fa 0.0112 the branch turns again towards increasing fi, now with two unstable real eigenvalues 
[Fig- Et c )]- With further increase in \i the L 2 norm of these unstable states continues to decrease as 
the solution begins to resemble more and more a 2d steady flowing film whose profile is modulated 
by the hydrophobic defects in the substrate below. 

Since each of the three SR states at /i = can itself be continued in ji we expect to find 
additional SR states when ji > 0. Of these SR2 and SR3 annihilate in a saddle-node bifurcation at 
/ig„ 3 w 0.0065 [Fig. [TJa)] leaving only the state SRi at larger fi. These states play no role in the 
depinning transition and the bifurcation at fj^ 3 is a regular saddle-node bifurcation. 

For larger ridge volumes (e.g., h = 1.3) the bifurcation diagram simplifies since the Hopf 
bifurcation is now absent. In this case the branch of stick-slip ridges (SSR) emerges directly from 
the saddle-node bifurcation in a sniper bifurcation, as confirmed by the fact that the inverse time- 
period now varies as (/i — /igni) 1 ^ 2 near me bifurcation. For a more detailed discussion of the 
transition between these two scenarios see Ref. [391. 




2. Transverse instability of a ridge 

A real 3d system will, however, show the behavior described in the previous section only if 
the spanwise domain size L is small. For large or indeed infinite L the ridge solutions (2d drops) 
are unstable, as already noted, with respect to spanwise perturbations. On a horizontal substrate 
(// = 0) the only relevant instability of this type is the Plateau-Rayleigh instability discussed in 
section HTlAI In this section we discuss analogous instabilities when /i > 0. 

To do so we follow the SRi branch in Fig.|7£a) as a function of fx together with the eigenvalues 
determining its linear stability properties with respect to spanwise perturbations of spatial period 
L. This procedure allows us to identify, as a function of \i, the critical spanwise domain size L c 
for the onset of a linear instability. The underlying linear stability analysis is performed using the 
Ansatz h(x, y) = h (x) + ah±(x) exp(ik y y) with k y = 2n/L. 

The resulting stability diagram in the (L, fi) plane is shown in Fig. [TOl Figure [TT] shows four 
sample perturbation profiles hi(x) along the stability boundary together with the associated ridge 
profiles h (x). These are discussed below. 

On the horizontal substrate (fi = 0) the SRi states are stable below the critical length L° cl . 
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FIG. 10: (color online) Stability diagram for steady ridge states on a heterogeneous incline showing the 
region of linear stability in the (L, /i) plane. Solid (dashed) lines indicate the presence of 3d (2d) instability. 
The remaining parameters are as in Fig. [7] 

When fi > these states become asymmetric with respect to the reflection x — > —x [Fig. [Ztb)] 
but remain stable for \x < nlf(L). For l? cl < L < L max the SRi states are stable or unstable 
depending on the value of fx. This is a consequence of the nonmonotonic dependence of L ci on fi 
shown in Fig. [10] for e = 0.3: the value of L cl first increases from L° cl to a maximum value L max 
at /x max before decreasing towards L min at fi^ni- F° r me present parameter values L min < L° cl and 
linear considerations alone allow one to distinguish four qualitatively different responses to the 
driving force fi, depending on the lateral system size L: 

(i) For L < L min = 22.4 the pinned ridge is linearly stable with respect to spanwise per- 
turbations in its entire range of existence, i.e., for all /i < /i^. The depinning behavior 
corresponds to the 2d case. 

(ii) For Lmin < L < LP cl the ridge is linearly stable when /i = but loses stability at a driving 
force fj%i < ulni- The instability has largest amplitude on the downstream side of the ridge 
[Fig.fTTTa)] indicating the onset of spanwise modulations that invade the wettability defect. 
As L approaches L min , fi^f approaches /i^. 

(iii) For L° cl < L < L max the SRi state is linearly unstable even when /i = 0. Increasing fi 
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X X 

FIG. 11: (color online) The marginal eigenfunctions h\{x) along the stability boundary in Fig. [11)1 (upper 
panels) and associated steady ridge profiles ho(x) (lower panels) at (a) fj, = 0.01157, L = 22.5; (b) 
H = 0.01019, L = 28.0; (c) n = 0.00147, L = 28.0; (d) ,u = 0.00656, L max = 30.1. 

stabilizes the ridge at a critical value /^f- The marginal eigenfunction [Fig. [TTT c)l peaks 
slightly on the upstream side of the ridge indicating stabilization with respect to the Plateau- 
Rayleigh instability. The ridge is then linearly stable up to the critical value nlf, where it 
loses stability as in case (ii), cf. Fig. fTTT a.b). 

(iv) For L > L max = 30.7 the ridge is linearly unstable for all ji. However, linear analysis is 
not able to tell whether the ridge will evolve into steady pinned drops or sliding drops. We 
expect that a critical ji exists below [above] which the former [latter] occurs. 

The different types (i)-(iv) of linear behavior inevitably result in different nonlinear behavior. 
For example, in cases (ii) and (iii) depinning occurs in an intrinsically 3d manner, as individual 
"fingers" extend across the defect. Thus the SRi state does not slide over the defect as a whole. 
In the following section we analyse the relation between drop and ridge solutions in the parameter 
regions (i) to (iv). We find that the intricate nonlinear behavior that results requires the introduction 
of several qualitatively different subregions in parameter space. These are discussed in detail in 
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the next Section before an overview of all stable solutions is given in Section ITlIDl 
C. Depinning of drops and ridges on an inclined substrate 
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FIG. 12: (color online) Bifurcation diagram showing the L 2 norm \ \Sh\\ of steady solutions SR2 and SR3 
as a function of the driving force fi together with the bifurcating transversally modulated ridges SD4. The 
latter branch is shown for the two different values L = 17.9 (red dotted) and L = 16.0 (blue dashed). The 
number of red +'s indicates the number of unstable eigenmodes when L = 17.9. Remaining parameters 
are as in Fig. |2j 

The analysis of the linear stability of the steady ridge solutions for different spanwise system 
sizes performed in the previous section indicates that the relation between drop and ridge solu- 
tions changes with lateral system size, i.e., with drop volume. The linear analysis allowed us to 
distinguish four parameter regimes. The curve of neutral stability w.r.t. harmonic spanwise pertur- 
bations in Fig. [10] indicates the loci of bifurcation points where modulated ridge solutions emerge 
from the spanwise-invariant ridges. In the following, we continue these modulated ridge solutions 
for several fixed values of L while changing fi. In addition, we continue the steady drop-like states 
present at jj, = (section ITU A[ Fig. [5]) towards larger ji. All solutions obtained in this way are pre- 
sented in a sequence of bifurcation diagrams together with the spanwise-invariant ridge solutions. 
These show branches of time-periodic solutions, i.e., drops and ridges that either oscillate or slide 
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from defect to defect, in addition to the steady solutions. 



1. Scenario (i): L < L m j n 

For L < L min the pinned ridge is linearly stable w.r.t. spanwise perturbations and depins via 
a Hopf bifurcation at /x^opf < f^sn- The marginally stable eigenfunction at the Hopf bifurcation 
remains translation-invariant in the spanwise direction. Consequently the depinning scenario in 
this case is identical to that already described for the 2d case in section IIIIB II and summarized 
in Fig. [7 
Refs. [139, 

become unstable to three-dimensional perturbations at some ji > /jl^ 4 , resulting in a sliding 3d 
state. 

The secondary branch SD 4 present on the horizontal substrate for L° 3 < L < LP c2 < L min 
(Fig. [5]) plays no role in any of the depinning scenarios discussed below. Continuation in fi shows 
that this solution terminates either on SR 2 or SR 3 , depending on the value of L (FigJT2"l). The SD 5 
solution likewise plays no significant role despite its presence in most of the scenarios discussed 
below (cf. Fi g s . [T51 and 1251 below) . 



After depinning the ridge undergoes stick- slip motion as described for 2d drops in 
, l^]. Depending on the value of L < L min this 2d time-dependent state may in turn 



2. Scenario (ii-a): L m i n < L < L^ nl 

For L min < L < lP cl the pinned ridge SRi is linearly stable w.r.t. spanwise perturbations at 
small fx but becomes unstable w.r.t. 3d perturbations at some nlf(L) < /i^. Thus the Hopf 
bifurcation that produces depinning may be preceded by a steady state spanwise instability. In this 
case the depinning that takes place at larger fi corresponds to an instability of a 3d steady state. 

We have already seen that when \i = and L < L° nl we have the three SR solutions shown 
in Fig. [7] For L® nl < L < LP cl we have in addition two 3d solutions, one ridge-like and the other 
drop-like. Each of these solutions extends smoothly to fi > with both L snl and L cl changing 
with fi, with L sn i(/i) < L° nl . Thus for L mill < L < L° nl the bifurcation diagram takes the form 
shown in Fig. [13] 

Here the branch of largest norm (thin line) corresponds to the steady spanwise-invariant ridges 
SRi. This solution loses stability with respect to 3d perturbations at /j^f w 0.01155. This value 
is slightly larger than /Jw, f ~ 0.01154 implying that the SRi state depins via 2d depinning. The 
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FIG. 13: (color online) Scenario (ii-a): (a) Bifurcation diagram for ridge and drop solutions on a substrate 
with hydrophobic line defects (wettability contrast e = 0.3) and spanwise system size L mm < L = 23 < 
Lg nl showing the L 2 norm \\Sh\ \ of steady solutions as a function of the driving force /x. A zoom is given 
in panel (b). Both steady spanwise invariant ridges (SR, thin red lines), and secondary drop solutions 
(SD2 and SD3, black lines) are included. Solid [dashed] lines indicate linearly stable [unstable] solutions. 
Downward [upward] pointing triangles indicate SSR [SSD] solutions. The former are taken from Fig. |7jb) 
and are stable with respect to 3d perturbations for fi < 0.015 only. Square symbols indicate solutions whose 
profiles are shown in Fig.|T4j The SD3 branch acquires stability at /^opf ~ 0.092 (off scale). The remaining 
parameters are as in Fig. [3] 

3d steady states produced at /xj?i bifurcate towards smaller /x before turning around towards larger 
/1 at a saddle-node bifurcation at /x^ « 0.005532 and are unstable throughout. In the follow- 
ing we use SD 2 to refer to the solutions above the saddle-node and SD 3 to refer to those below. 
As shown in Fig. [14] the appearance of the 3d solutions changes substantially along the branch. 
Near /j^f [Fig. [T4T a)l the solution is ridge-like. At the saddle-node at /x^i it becomes drop-like 
[Fig.[T4Tb)l while at large /x the solution becomes a steady (unstable) streamwise-modulated rivulet 
[Fig. (He)]. 

As already mentioned the 2d SRi state loses stability at /x = /^o P f with respect to 2d oscillations 
(see section IIII B 2[ Fig. [7]) and these must undergo a global bifurcation resulting in depinning 
(section ITVl). The resulting SSR state remains stable for a range of values of /x but loses stability 



23 




(a) (b) (c) 

FIG. 14: (color online) Scenario (ii-a): Unstable steady states at locations indicated by open squares 
in Fig. [13] in terms of contours of constant h(x,y). (a) Drop solution SD2 at ll = 0.0105 resembling a 
modulated ridge, (b) SD2 solution very close to the saddle-node at fi^i ~ 0.005532. (c) Static rivulet SD 3 
at fi = 0.0253. Parameters are as in Fig. |T3l 

at /i PS 0.02 to spanwise perturbations resulting in a hysteretic transition to stable stick-slip drops 
(SSD). The latter appear to undergo a saddle-node bifurcation near ll = 0.0123 and do not reach 
any of the steady states at lower ll. 

The unstable steady rivulet states on the SD 3 branch [Fig. [LU c)! stabilize through a Hopf 
bifurcation at a much larger /i, //w>f ~ 0.092. If a time simulation is done for /1 < fJ^ p{ (but 
above /i ~ 0.01) with the unstable steady rivulet solution as the initial condition the rivulet decays 
into sliding drops corresponding to a state on the SSD branch. Thus the SSD branch terminates in 
a Hopf bifurcation on the SD 3 branch (cf. Fig. [251 below). 

3. Scenario (ii-b): L® nl < L < L° 

For L > L° nl two 3d solutions are present (in addition to the previously studied 2d SR solu- 
tions) even at ll = 0. Another two appear at L = L° n2 slightly above L° nl (cf. Fig. [5]). The resulting 
bifurcation diagram is shown in Fig. [15] As in scenario (ii-a), the pinned 2d ridge is linearly stable 
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FIG. 15: (color online) Scenario (ii-b): Bifurcation diagram for ridge and drop solutions on a substrate with 
hydrophobic line defects (wettability contrast e = 0.3) and spanwise system size L = 26 showing the I? 
norm of steady solutions as a function of the driving force \i. Both steady spanwise invariant ridges 
(SRi, thin red lines), and secondary drop solutions (SDi, SD2, SD3 and SD5, black lines) are included. 
Solid [dashed] lines indicate linearly stable [unstable] solutions. Downward [upward] pointing triangles 
indicate SSR [SSD] solutions. Solution profiles at locations indicated by open squares are shown in Fig. [18] 
The remaining parameters are as in Fig. [3] 

at small \x but becomes unstable w.r.t. 3d perturbations at ii z &{V) < /i^Lf < //gni- Continuation 
of the branch of SD 2 solutions when L = 26 shows a branch that bifurcates subcritically from the 
SRi branch at fj^f ~ 0.0111 (marked by a cross in Fig. [T31) and continues all the way to ji = 0. 
The two solutions with the largest norm at fi = correspond to stable and unstable steady drop 
(SDi) solutions, connected via a saddle-node bifurcation at /i^ ~ 0.002019. Such solutions are 
present whenever L > L° nl . Time simulations at [i slightly above {1^2 w i m a stable steady drop as 
initial condition show that the saddle-node bifurcation is not a sniper bifurcation: the steady drop 
does not start to slide but instead stretches in the y-direction and converges to the stable 2d SRi 
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FIG. 16: (color online) Temporal period T of stick-slipping drop (SSD) solutions as a function of the driving 
force /i for L = 26 and L = 27 (see legend). In each case the first six points were used for a fit of the form 
T = — aoln(/x — a\) (dashed lines). The fit parameters are ao = 88.3, a\ = 0.00939 for L = 26, and 
a = 112.85, a x = 0.00810 for L = 27. 



state. 



The transition from (ii-a) to (ii-b) is now clear: with increasing L the saddle-node bifurcation 



between the branches SD 2 and SD 3 [at /i^, cf. case (ii-a)] moves towards yt, — 0. At L — L° 



sn2 



it touches the line /i = and two new solutions appear at [i — corresponding to the SD 2 and 
SD 3 branches in Fig. [5]). In contrast, the behavior of the sliding solutions remains unclear. The 
2d oscillations and stick-slip states (Fig. [7]) continue to exist but are now unstable with respect 
to 3d perturbations, and time- simulations result in stick-slip drops (SSD). Since this branch does 
not terminate in a local bifurcation we conjecture that it terminates in a global bifurcation on the 
(unstable) SD 2 branch. The discussion in section [IV] supports this interpretation as does Ref. |}86fl - 
In Fig. [16] we show logarithmic fits to the SSD period T, i.e., the time it takes for a drop to slide 
to the next downstream defect following depinning, as a function of /i for L = 26 (cf. Fig. [T51) 



and L = 27 (cf. Fig. [T9l) . The points near the critical parameter value fj, — ax « ^ d are given 
the highest weight. A square-root power law of the type expected near a sniper bifurcation is not 
compatible with the data while the logarithmic fit shown in the figure provides strong evidence for 
depinning via a global bifurcation for both L = 26 and L = 27. We caution, however, that fits of 
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FIG. 17: Real part of the leading eigenvalues of the SD2 and SD3 branches for (a) L = 26 (Fig. |T31) and (b) 
L = 27 (Fig. O. The dotted line indicates the approximate location of the global bifurcation, i.e., fj, = a\ 
(Fig. [16]). As /i increases the SD3 branch stabilizes via a Hopf bifurcation at /^ pf ~ 0.075 (L = 26) and 
at fj,^ f « 0.07 (L = 27), respectively. 

this type cannot exclude the presence of a fold in the SSD branch prior to depinning (fi — > /i^ d , 
T — > 00) because of the absence of very high period simulations. On the other hand, the figure 
clearly differentiates between the global bifurcations present for L = 26 and L = 27 and sniper 
bifurcations present for L = 27.8 (cf. Fig. [22 below) and larger L, for which a square-root fit 
works very well. 

The logarithmic fits provide information about the leading eigenvalues of the saddle-type equi- 
libria responsible for the presence of the global bifurcation suggested by the fits. We assume that 
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(a) (b) (c) (d) 

FIG. 18: (color online) Scenario (ii-b): Steady states at locations indicated by open squares in Fig. [T5l in 
terms of contours of constant h(x, y). Drops at (a) the saddle-node /i^ 2 « 0.0020, and on the SD 2 branch 
at (c) /j, = and (d) fi = 0.0070. The modulated ridge shown in (b) is also at \i = 0. 

the two leading eigenvalues of the equilibrium at fj, — a\ (i.e., the eigenvalues closest to zero) 
are real, with one positive eigenvalue \ u > and one negative eigenvalue \ s < 0. We assume 
that all the remaining eigenvalues, whether real or complex, are stable with modulus larger than 
|A 8 |. Under these circumstances standard theory shows that near fj, — a%, the period follows the 
asymptotic behavior T — j- ln(/i — a%) provided X u + X s < 0, and T ps j- ln(/i — a x ) provided 
A u + A s > 0. In the former case the periodic orbit is stable near fi = a,\, in the latter it is unstable. 
However, in both cases the growth of the period is asymptotically monotonic with the distance 
from n = a±. 

We have computed the leading eigenvalues along the SD 2 and SD 3 branches for the two val- 
ues of L. Both branches are of saddle-type near the global bifurcation and hence candidates for 
involvement in the global bifurcation. Figure [FT] shows that for (a) L = 26 the two leading eigen- 
values of SD 2 at n = a x « 0.00939 are A u = 0.0140 and A s = -0.00396, while for (b) L = 27 
the leading unstable eigenvalue at fi = a\ ps 0.00810 is \ u = 0.0101 and the leading stable eigen- 
value is complex, with real part ReA s = —0.00420. Thus the theory summarized above applies 
in the first case only; the complex eigenvalues in the second case suggest an oscillatory approach 
to the global bifurcation of Shil'nikov type that is not observed. In the first case the fact that 
A u + A s > indicates that the appropriate prefactor in the scaling of the asymptotic period is 
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ao = 1/| A s | ~ 252.5. This value does not compare well with the fit (a) in Fig. [[6] As in the 2d 
case we surmise that this discrepancy is due to the fact that X u + A s > implies that if SD 2 is 
responsible for the global bifurcation the nearby periodic orbits will necessarily be unstable. Thus 
the branch of periodic states must in fact overshoot the global bifurcation, and double back, i.e., 
the branch must undergo a saddle-node bifurcation and lose stability prior to the global bifurcation 
much as in the 2d case discussed in Fig. U\ However, in contrast to Fig. U\ in the 3d case we have 
been unable to detect any evidence for a saddle-node of the SSD states in this region. Likewise in 
the second case the fact that X u + 2ReA s > implies that all long period SSD states will also be 
unstable. 

Close to ji = fig d , but before the asymptotic regime just described sets in, the periodic orbits 
are, however, stable, and in this case the period should be approximated by T « — y- ln(/i — ai), 
i.e., we expect that in this regime the correct value of ao is oo ~ 71.4 (L = 26). This result is to 
be compared with the fit a = 88.3 (Fig. [16]). Of course, because of the overshoot, the fitted value 
of ai, ai »j 0.00939, is smaller than the true value ji = [i 3 g d and this is so in case (b) as well. Here 
the fit yields a « 112.85 (L = 27) which again compares well with j- = 101.3. 

4. Scenario (iii-a): l? c < L < L\ 

As L increases beyond L° c two related changes occur simultaneously at [i = 0: the 2d SRi 
state becomes linearly unstable and the subcritical part of the SDi branch is no longer present 
(cf. Fig. [3]). As a result when \l increases the SRi acquires stability with respect to the Plateau- 
Rayleigh instability at /i™ « 0.000286 (for L = 27) before losing it again at /i 3 c f w 0.0107783 
(Fig. [TO]). The resulting bifurcation diagram obtained using continuation together with time- 
stepping computations is shown in Fig.[[9] Selected steady profiles are shown in Fig. [201 

Although the linear stability properties of the SRi state do not change when L increases towards 
L max , the bifurcation diagram changes its appearance twice owing to codimension two bifurcations 
that take place at L = L\ and L = L 2 (see below). We treat the resulting L regimes in cases (iii-a), 
(iii-b) and (iii-c), respectively. 

The main difference between case (iii-a) and the scenario (ii-b) already discussed is found in 
the behavior of the unstable part of the SDx branch. In (iii-a) this branch does not extend to 
fj, — but terminates instead on the SRi branch in a subcritical bifurcation at /i^ • With increasing 
amplitude the unstable part of the SDx branch turns around in a 'true' saddle-node bifurcation at 
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FIG. 19: (color online) Scenario (iii-a): Bifurcation diagram for steady ridges (SR, thin red lines), steady 
drops (SD, black lines) on a substrate with hydrophobic line defects (wettability contrast e = 0.3) and 
spanwise system size L = 27 showing the L 2 norm \ \5h\ \ of steady solutions as a function of the driving 
force fi. Solid (dashed) lines denote stable (unstable) solutions. Downward [upward] pointing triangles 
indicate SSR [SSD] solutions. Solution profiles at locations indicated by open squares are shown in Fig. [20] 
The remaining parameters are the same as in Fig. [3] 

A*sn2 ~ 0.0037015 and acquires stability. Time evolution at \i slightly above the saddle-node at /ij^ 2 
starting from a stable drop at lower /i converges to the stable SRi solution. Examples of profiles 
on the SDx branch are shown in Fig. [201 The bifurcation at ^ also generates subcritical 3d steady 
states. As in scenario (ii-b) these extend all the way to fi = 0, albeit with a slightly higher and 
more pointed maximum norm. The norm at small /i is also larger and is now comparable to the 
amplitude of the SRi states. For /i > trme evolution results in a family of stick-slip drops 
(SSD). 
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FIG. 20: (color online) Scenario (iii-a): Steady states at locations indicated by open squares in Fig. [19] in 
terms of contours of constant h(x, y), ordered by decreasing norm, (a) \i = 0.00199, (b) \i = 0.00203, (c) 
fx = 0.00188 and (d) \i = 0.00208. 

5. Scenario (iii-b): L\ < L < L 2 



At L = L\ = 27.59 a codimension two bifurcation takes place and the bifurcation diagram 
changes from Fig.[l9]to that in Fig.[2Ha). The figure shows that no qualitative change in behavior 
occurs at zero or small fj,. The SRi branch is unstable at \i = and stabilizes at fj^ ~ 0.0010798 
(for L = 27.8). This bifurcation continues to be subcritical and produces a branch SDi of unstable 
3d states. The SDi states annihilate in a 'true' saddle-node bifurcation at /i^ 2 w 0.0058161. 
However, the 3d solutions that emerge subcritically at fj^f fa 0.0104789 and extend to fi = 
undergo a major change near maximum norm. At L = L\ a cusp appears and for L > L\ 



the branch develops a loop, with two additional saddle-node bifurcations, at \\f s 



3d 
sn3 



0.0069872 



and /i^ iper = /4n4 ~ 0.0072342 (for L = 27.8), as visible in the zoom in Fig. EHb). At the 
same time the branch acquires a linearly stable segment that extends between the two saddle-node 
bifurcations. A sample stable solution is shown in Fig. [22} Moreover, while the left-hand saddle- 
node at //^g is a 'true' saddle-node, that at /J^ pei = fJ^ 4 now corresponds to a sniper bifurcation, 
with a branch of stick-slip drops (SSD) emerging from /i^pcr 

With increasing L the topology of the diagram remains unchanged until a second critical value, 
L = L 2 , as discussed next. 
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FIG. 21: (color online) Scenario (iii-b): (a) Bifurcation diagram for steady ridges (SR, thin red lines), and 
steady drops (SD, black lines) pinned by a hydrophobic line defect of strength e = 0.3 for L = 27.8 
showing the L 2 norm \\Sh\ \ of steady solutions as a function of the driving force fi. Solid (dashed) lines 
denote stable (unstable) solutions. The panel (b) shows a zoom of the cusp-like feature on the SD2 branch. 
The SD3 branch becomes stable at \i = 0.064 (not shown). Downward [upward] pointing triangles indicate 
SSR [SSD] solutions. Solution profiles at locations indicated by open squares are shown in Fig. |22] The 
remaining parameters are the same as in Fig. [3] 

6. Scenario (iii-c): L2 < L < L max 



Increasing L beyond L 2 = 27, 87 another codimension two bifurcation takes place and the bi- 
furcation diagram changes from Fig. 1217 a) to that in Fig. [23} Once again the branches at \i = and 
their stability properties remain unchanged, as do the properties of the SRi states. However, the 
nonlinear behavior at finite fj, changes dramatically with the 3d branches in Fig. l2TT a) reconnecting 
in a new way. This reconnection occurs as a result of the fusion/annihilation of the saddle-node at 
yu^ 2 on the SD 1 branch and the newly created saddle-node at ^3 on the branch of 3d states SD 2 
via a "necking" bifurcation. As a result the branches reorder into a simpler bifurcation pattern 
[Fig. l237 a)1. The stable part of the SDi branch now annihilates with the left part of the unsta- 



ble SD 2 branch in a saddle-node bifurcation at /isn ipcr = AW ~ 0.0069055 that is still a sniper 



bifurcation [in Fig. !23T a)1. As a result of the various transitions occurring with increasing L the 
depinning transition has shifted from the 2d SRi states [|39ll to the 3d SDi states 114 ill , i.e., the 
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FIG. 22: (color online) Scenario (iii-b): Steady states at locations indicated by open squares in Fig. [2T1 in 
terms of contours of constant h(x,y). (a) \i = 0.00082 on the drop branch SDi, (b) on the drop branch 
SDi at = 0.00110, (c) at the saddle-node /i^ 2 = 0.00582 on SDi, and (d) at the sniper bifurcation at 
/i^ = 0.00724 on SD 2 . 

behavior expected of large drops in large domains. 

Figure [23Tb ) shows the eigenvalues along the stable and unstable parts of the SDi branch and 
reveals that the leading eigenvalue of the stable drop solution is complex close to the saddle-node 
at /^sniper- m m i s regime the stable SDi states relax in an oscillatory fashion when perturbed. In 
contrast, the leading unstable eigenvalue along the unstable SDx branch is real, implying mono- 
tonic growth. This eigenvalue is created by a collision of the complex eigenvalues on the negative 
real axis prior to the saddle-node bifurcation creating a pair of real eigenvalues, one of which 
crosses into the positive half-plane at ^ ipei . 

At the same time the branch of unstable 3d states that emerges at fj^ ps 0.0025782 continues 
towards and ends at fj^f « 0.010065. As L is increased further, an d approach each other 
until they annihilate at L max (Fig.flOl). Consequently the branch SD 2 of unstable 3d states becomes 
shorter and shorter, along with the interval of stable SRi states, and when it vanishes so does the 
SRi stability interval. However, for the parameter values of Fig. [23^ a) we find coexistence between 
stable SRi states and either stable SDx states or, beyond the sniper bifurcation at l^ ipev = ft 3 ^ 



(numerically stable) stick-slip drop motion (SSD). 
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FIG. 23: (color online) Scenario (iii-c): Bifurcation diagram for steady ridges (SR, thin red lines), and 
steady drops (SD, black lines) on a substrate with hydrophobic line defects (wettability contrast e = 0.3) 
and spanwise system size L = 29 showing the L 2 norm \\Sh\\ of steady solutions as a function of the 
driving force fj,. Solid (dashed) lines denote stable (unstable) solutions. Downward [upward] pointing 
triangles indicate SSR [SSD] solutions. Solution profiles at locations indicated by open squares are shown 
in Fig. [24] The remaining parameters are as in Fig. [3] The SD3 branch becomes stable at \i = 0.062 (off 
scale), (b) Real part of the leading eigenvalues for the stable (thin black line) and unstable part (heavy red 
line) of SDi drop branch. The black dashed line corresponds to a complex eigenvalue. 



7. Scenario (iv): L max < L 

When L increases beyond L max the linear stability problem indicates (Fig. [TO]) that the 
spanwise-invariant ridge SRi is linearly unstable for all /i with a real positive leading eigenvalue. 
Thus the two pitchforks responsible for the unstable modulated ridges annihilate as L increases 
and the modulated ridge states disappear (see Fig. [25]). As a result the 3d drop states SDx now 
correspond to the drops studied for large domains in J^j]]. At these values of L a stable and an 
unstable drop of this type annihilate in a sniper bifurcation at /i^ 4 w 0.0069. The emerging branch 
of stick-slipping drops (SSD) gradually approaches the branch of SD 3 states, connecting to it in 
a Hopf bifurcation at /x^opf ~ 0.051. Since at large /i the SD 3 states resemble modulated rivulets 
[cf. Fig. [14] (c)] close to the Hopf bifurcation the SSD states correspond to surface waves on a 
steady rivulet. 
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FIG. 24: (color online) Scenario (iii-c): Steady states at locations indicated by open squares in Fig. [23] in 
terms of contours of constant h(x, y), ordered by decreasing norm, (a) \i = 0.00419, (b) \x = 0.00391 and 
(c) \i = 0.00401. The remaining parameters are as in Fig. [3] 

D. Stability in the (L, fi) plane 



The bifurcation analysis in the preceding sections has identified several types of stable steady 
(ridge, drop and rivulet) and time-periodic (stick-slipping ridge and stick-slipping drop) states. 
Figure [26] displays the stability regions of these states in the (L, fi) parameter plane. Most of the 
boundaries of these regions are computed by numerically tracking the bifurcations that lead to 
loss of stability. For example, since the drop state loses stability via a saddle-node or a sniper 
bifurcation, one of the stability boundaries requires the tracking of a zero eigenvalue of the linear 
problem obtained by linearizing Eq. (OQ) about the drop state. This is accomplished using a steady- 
state continuation algorithm similar to that described in [54]. Likewise, since the rivulet state is 
stabilized via a Hopf bifurcation numerical tracking of the Hopf bifurcation yields the stability 
boundary of the rivulet state. 

The region of stable steady ridge solutions is discussed in section HIlB 21 Note that it has some 
overlap with the region of stable steady drop states. This indicates bistability and implies that, in 
the intersection region the observed state depends on initial conditions. The region of stable steady 
drop solutions is delimited by two curves which end in a cusp (upper green triangle) corresponding 
to the codimension two bifurcation at (Li, //i) (see section UlIC 51) . The leftmost boundary of the 
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FIG. 25: (color online) Scenario (iv): Bifurcation diagram for steady ridges (SR, thin red lines), and steady 
drops (SD, black lines) on a substrate with hydrophobic line defects (wettability contrast e = 0.3) and 
spanwise system size L = 32 showing the L 2 norm \ \5h\ \ of steady solutions as a function of the driving 
force Solid (dashed) lines denote stable (unstable) solutions. Downward [upward] pointing triangles 
indicate SSR [SSD] solutions. The remaining parameters are as in Fig. [3] The SD3 branch becomes stable 
at < pf « 0.051. 

SD region (thick blue dotted line in Fig. [26]) starts from the point (L° snl , \x = 0) and terminates at 
the cusp. The turning point at (L 2 , /i 2 ) (lower green triangle) corresponds to the codimension two 
bifurcation identified in section UlI C 6[ This dotted boundary is the locus of 'true' saddle-nodes, 
i.e., a drop destabilized at this boundary does not start to slide but instead stretches and converges 
to a stable ridge solution (see section llJIC 31) . In contrast, the upper boundary of the SD region (red 
dashed line in Fig.l26l) corresponds to sniper bifurcations, and loss of stability along this boundary 
results in periodic stick-slip motion. This scenario occurs even if the sniper boundary intersects 
the stable ridge region. The resulting stick-slip motion is stable as shown by time integration, and 
corroborated in the next section. 
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FIG. 26: (color online) Stability diagram for steady (SR) and stick-slipping ridge (SSR), steady (SD) and 
stick-slipping drop (SSD) and rivulet states on a heterogeneous incline showing the region of stability in the 
(L, fi) plane for e = 0.3, h = 1.2, and L x = 40. Black solid (dashed) lines indicate the transversal (sniper) 
instability of the steady ridges as in Fig. [10] The thick dotted blue (dashed red) lines are the loci of saddle- 
node (sniper) bifurcations of the drop states. The thin dotted blue line indicates the hypothetical border of 
the region of stable stick-slipping drops (the three "+" symbols result from our calculations). Finally, the 
dot-dashed green line indicates the loci of the Hopf bifurcation where the steady 3d rivulets become stable 
when n increases. 



We are not, however, able to compute in detail the region of stable stick-slip motions for L < L\ 
since this requires the knowledge of the loci of global homoclinic bifurcations, a difficult task for 
a problem of this complexity. Instead we indicate the hypothetical boundary of the region of stick- 
slip drop motion by a thin blue dotted line, guided by the orange crosses. These indicate the 'end' 
of the branch of stick- slipping drops as found in time simulations for selected L. Note, however, 
that the transition at the thin blue dotted line involves some hysteresis, at least outside the steady- 
ridge region, where simulations indicate some overlap between regions of stable stick-slip ridges 
and stick-slip drops. The extent of this hysteresis region remains unknown. 

At large ji one encounters a region of stable rivulets. This region does not intersect the other 
stability regions (note the 'broken' y-axis in Fig. [26b. The critical driving jj hop f(L) corresponding 
to the stability boundary increases rapidly with decreasing L. 
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IV. INTERPRETATION 



In this section we provide a bifurcation theory interpretation of the transition sequences de- 
scribed in the preceding section. We focus on the transitions observed for L = 29, L = 27 and 
L = 23 as the forcing fi varies. These transitions involve first and foremost the saddle-node bi- 
furcation^) on the steady ridge branch SRi (hereafter SR). In addition to this bifurcation we have 
identified two nearby bifurcations, a Hopf bifurcation that leads to time-periodic oscillations of 
the ridge state, and a steady state bifurcation that breaks the invariance of the ridge state in the 
spanwise direction, leading to steady drop-like states (SD). In addition, we have found two time- 
dependent solution types, sliding stick-slip ridges (SSR) and sliding stick-slip drop-like states 
(SSD). 

Our computations reveal that both the Hopf and symmetry-breaking steady state bifurcations 
occur close to the upper saddle-node bifurcation, suggesting that the different bifurcation scenarios 
may be understood by examining the interaction between these three bifurcations, i.e., a codimen- 
sion three bifurcation. The inclusion of the nearby lower saddle-node bifurcation would require a 
study of codimension four in which a hysteresis bifurcation interacts with both a Hopf bifurcation 
and a steady-state bifurcation. Bifurcations of this complexity have not been studied in full detail, 
eg. [1870. Moreover, the standard approach using normal form theory describes solutions in R 3 , the 
dimension of the (extended) center manifold for a mode interaction of this type. This reduction in 
dimension is a consequence of normal form symmetry which decouples the phase associated with 
the Hopf bifurcation from its amplitude. 



89] 



The three-dimensional phase space R 3 allows of course complex (chaotic) dynamics II88L 
but is unable to accomodate the sniper bifurcations that lead to sliding states. For this purpose 
the problem must be posed in R 2 xS x , i.e., with periodic boundary conditions inherited from the 
periodicity of the problem in the downstream direction. This type of phase space allows two types 
of periodic orbits that oscillate about an equilibrium (librations) and those for which the periodic 
variable increases by 2ir each period (rotations). These solution types can be visualized in terms 
of the standard pendulum, with librations corresponding to standard pendulum oscillations and 
rotations corresponding to a whirling pendulum. 

Although no rigorous normal form theory exists for this case (the rotations rely on global 
properties of the phase space) model problems with the required properties are readily constructed 
(eg. The most relevant such model was constructed by Krauskopf and Oldeman [90] who 
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investigated the properties of the following system of equations: 



P 



v\p — ap sin <p — p 3 , 



(6) 



ip = v 2 + sp 2 + 2 cos <p + cp 4 . 



(7) 



The model can be viewed in two ways, as an interaction between a saddle-node bifurcation in a 
periodic orbit (sniper) and a Hopf bifurcation, or an interaction between a sniper and a symmetry- 
breaking pitchfork bifurcation. In the former the variable p > corresponds to the amplitude of 
the periodic oscillations (the temporal phase decoupling in "normal form") while in the latter p 
corresponds to the amplitude of the state created in the symmetry-breaking pitchfork bifurcation. 
Both interpretations rely on the symmetry p — > — p of the model and are relevant to the depinning 
problem. 

We examine first the depinning problem in which the steady ridge state (SR) does not undergo 
a Hopf bifurcation. In this case the system ©-([7J) describes the dynamics near the codimension 
two point (per, Lct) at which the symmetry-breaking pitchfork of the SR state coincides with 
the saddle-node bifurcation at p = p 2 ^. Thus (u 1: u 2 ) represent linear combinations of p — p, C T 
and L — Lct- In this case the variable p > represents the amplitude of the translation invariance 
breaking instability creating the SD state, the spanwise phase of the mode having decoupled be- 
cause of translation invariance of the SR state in the spanwise direction. The variable <p plays the 
role of the downstream coordinate. Thus equilibria with p = represent SR states while periodic 
solutions with p = correspond to sliding ridges (SSR). Equilibria with p ^ correspond to 
stationary drops (SD) while periodic solutions with p ^ are of two types, oscillations about an 
SD state (librations) and sliding drops (rotations, SSD). 

In Fig. [27] we reproduce Figure 7 of ||90|l . The central part of the figure shows the (ui, z/ 2 ) 
parameter plane; in the depinning problem these parameters correspond to different combinations 
of the parameters p and L. The (z/ 1; z/ 2 ) plane is split into a number of distinct regions, labeled by 
integers, with distinct phase space dynamics. The phase portraits characteristic of each region are 
shown along the periphery of the diagram, with the variable x shown horizontally and <p vertically. 
The lines So in the (ui, z/ 2 ) plane correspond to saddle-node bifurcations with p = 0, i.e., to snipers 
on the SR state, and are located at v 2 = ±2. The ellipse P corresponds to pitchfork bifurcations 
of equilibria with p = 0, i.e., to steady bifurcations from SR to SD. The line S corresponds to 
saddle-nodes of SD states (the point labeled BT is a Takens-Bogdanov point with a double zero 
eigenvalue). Finally the curves labeled Z indicate the transition from SSR states to SSD states, 
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FIG. 27: The (1/1,1/2) parameter plane for the system (O-© with a = 0.5, c = 0, s = — 1 (case B of 
|90|]) showing the codimension one lines So, P, S and Z discussed in the text. The Takens-Bogdanov 
point is labeled BT. The remaining curves correspond to global bifurcations. The phase portraits around the 
periphery show solutions characteristic of each of the numbered regions in the (p, cf>) plane. Reproduced 



with permission from [90]. 



i.e., transitions from a rotation with p = to a rotation with p > 0. The curves S , P and Z meet 
at points labeled B ± . The additional lines all correspond to different types of global bifurcations 
detailed in B2QD. 
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We are now in a position to interpret the dynamics observed in this paper in terms of the 
truncated "normal form" ©-CO). Specifically, the pinning transition observed as p decreases, for 
example, for L = 32 (Fig. [251), corresponds to traversing from region 13 (unstable SSR, stable 
SSD) to region 12 (a pair of unstable SR and a stable SSD) to region 3 (a pair of unstable SR and 
a pair of SD, one stable, one unstable). The transition from region 12 to region 3 is thus a standard 
sniper, except that the equilibria created in the transition are of SD type. 

When L = 27 the depinning transition summarized in Fig. [19] is consistent with a slightly 
different cut through the (u±, u 2 ) plane. This time we pass from region 13 to region 12 followed 
by regions 29 (a pair of SR states, one of which is stable and the other unstable, an unstable SD 
state created in a subcritical pitchfork bifurcation, and a stable SSD state) and region 22 (a pair 
of SR states, one of which is stable and the other unstable, and an unstable SD state). Note that 
the disappearance of the SSD state is the result of a standard homoclinic bifurcation, a collision of 



a periodic orbit with a single saddle-type equilibrium, here an SD state, as suggested by Fig. | 
Ref.H shows that these transitions are robust for a > 0, c = 0, s < [Eqs. ©-CO)] and represent 



the simplest transitions from among the possibilities shown in Fig. [27] and related figures valid in 
other parameter regimes. 

We next consider the transition observed for L = 23 (Fig. [13]). We have not been able to 
interpret this transition within the two-dimensional system ©-CO). The main obstacle for this 
interpretation is the coexistence of a stable SSR state with a pair of unstable SR states, a possibility 
that is not permitted within Eqs. ©-CO) when p — 0. Indeed, Fig. [13] indicates that in this case 
the SR state depins via a Hopf bifurcation, a possibility that can only be captured by a higher 
codimension normal form. On the other hand the observed transition from SSD to SSR resembles 
that along the line Z between regions 13 and 14 in Fig. [27] except for the fact that in Fig. [T37 b) 
this transition is apparently subcritical. 

Since the "normal form" ©-CO) wa s constructed to describe the interaction between a Hopf 
and a sniper bifurcation we can also examine its applicability in the two-dimensional depinning 
problem, i.e., for spanwise invariant ridge states. Figure[2b) shows the transition observed when 
three-dimensional instabilities are suppressed. The same transition forms part of the behavior in 
Fig. [13] despite the presence of 3d instability. To use Fig. [27] we need to recall that fixed points 
of Eqs. ©-CO) with p > now correspond to constant amplitude oscillations (OR) about an SR 
state. Thus the standard sniper bifurcation corresponds to a transition from region 14 to region 15. 
But when a Hopf bifurcation occurs on SR prior to the saddle-node at /i^ the sniper bifurcation 
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is replaced by a new depinning transition as described in §111 B I The SSR state now involves 
oscillations about the SR state and so is represented in Fig. [27] by a rotation at finite p. As seen 
from the behavior of the period in the inset in Fig.[7£b) this state disappears in a global bifurcation. 
We expect that this global bifurcation involves the unstable oscillations (librations) created at 
/i£o pf . In Fig. [27] this sequence of transitions is described by traversing from region 13 (stable 
SSR state with superposed oscillations) into region 12 (a pair of unstable SR states together with 
a stable SSR state with superposed oscillations) and then into region 29 (a pair of SR states, one 
of which is stable and the other unstable, together with an unstable libration and a stable SSR 
state with superposed oscillations). The pinning transition takes place on crossing from region 29 
into region 22 (a pair of SR states, one of which is stable and the other unstable, together with an 
unstable libration) and so involves the formation of a homoclinic connection to a periodic orbit. 
Note that in this scenario the SSR state in Fig.[7£b) is a two-frequency state while the libration that 
is left has a single finite frequency. 

Our calculations indicate that the SSR state is in fact a single-frequency state. As a result we 
favor a different explanation for the behavior shown in Fig. Hb). This behavior involves a global 
bifurcation between a single frequency rotation (SSR) and a single-frequency oscillation (OR) 
about an SR state, a transition that is organized by a Takens-Bogdanov bifurcation with periodic 
reinjection. At the Taken-Bogdanov point the Hopf bifurcation coincides with the saddle-node 
and the Hopf frequency vanishes. Consequently the "normal form" ©-CO) no longer applies. The 
Takens-Bogdanov point can be identified by shifting the parameter e appropriately and recomput- 
ing the leading eigenvalues along the SR branch to find the point (/itb, £tb) at which the zero 
eigenvalue has double multiplicity. Although we have not performed such a calculation we ex- 
pect this point to lie close to our chosen parameter values, implying that this bifurcation is likely 
relevant to the observed dynamics. The normal form for this bifurcation is also two-dimensional 
and we conjecture that in the presence of reinjection it exhibits the necessary gluing bifurcation. 
Similar bifurcations are known to occur in spontaneous parity breaking in inhomogeneous systems 
||9l|l . However, a study of this problem along the lines of [90(] is beyond the scope of this paper. 

The virtue of relating our results to something like normal form analysis is that the theory is 
able to identify "all" types of dynamical behavior that may be present in the vicinity of the origin 
in the (y\, u 2 ) plane [resp., (/i — p T B, e — ctb)L thereby indicating the possible dynamics that may 
be present in the depinning problem for other parameter values. On the other hand the restriction 
to a low-dimensional center manifold precludes the presence of complex dynamics. Such behavior 
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is typically associated with the presence of global bifurcations such as those taking place along 
the curves L, W and W in Fig. [27] In addition, the theory cannot account for bifurcations that 
occur "far" from the codimension-two point. Thus one does not expect to be able to explain all 
aspects of the computed bifurcation diagrams using low codimension normal form theory, and this 
is indeed the case. 



V. CONCLUSIONS 

In this paper we have studied the process of depinning of driven liquid drops and ridges on 
heterogeneous substrates and related transitions between drop and ridge states. In the absence 
of driving, e.g., on a horizontal substrate, the latter transition is related to the Plateau-Rayleigh 
instability of a liquid ridge. To explore the three-dimensional dynamics in the presence of driving, 
we have formulated the depinning process as a bifurcation problem, and focused on a generic 
problem of this type. 

We adopted a simple model problem in which the spatial heterogeneity corresponds to a modu- 
lation of the short range polar contribution to the disjoining pressure and takes the form of parallel 
hydrophobic stripes with a well-defined spatial period. The stripes are oriented transverse to the 
driving direction, and may block the passage of ridges and drops, producing strongly asymmetric 
ridges prior to depinning. In general, pinning might arise from spatially varying wetting properties 
as here, or from heterogeneous topography, temperature or electric fields. Possible driving forces 
include gravitational or centrifugal forces, and/or gradients of wettability, temperature or electric 
fields. In the present work we have used gravity as an example, i.e., we have studied drops and 
ridges on an incline. As a result our approach is limited to drops not much thicker than the wetting 
layer assumed to be present in our formulation. 

The two-dimensional version of this problem was formulated and studied in Ref. I39L where 
two possible depinning mechanisms were identified: depinning via a sniper bifurcation and depin- 
ning via a Hopf bifurcation followed by a global bifurcation that is required for translation of the 
ridge from one defect to the next. In the present paper we have been interested in like behavior 
for fully three-dimensional states that we variously refer to as modulated ridges or drops depend- 



ing on their appearance. Preliminary work on the 3d problem |41|] using the numerical procedure 
described in d4I1 indicated that for relatively large drops the 3d problem behaves qualitatively 
like the 2d problem but the detailed depinning mechanism was not studied owing to the long time 
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scales involved in the transition region. In this paper we have therefore focused on the details of 
the depinning transition in the 3d case. 

In three dimensions a transverse liquid ridge has additional avenues open to it since the whole 
ridge no longer needs to depin simultaneously but can instead send out protruberances through 
which the fluid flows over the heterogeneity. In this paper we have identified seven distinct regimes 
depending on the spanwise spatial period L of the computational domain. This spanwise length, if 
sufficiently small, can suppress vestiges of the Rayleigh-Plateau instability that leads to spanwise 
modulation of a ridge even on a horizontal substrate. However, on an inclined substrate the choice 
of L may not only allow a 3d instability to proceed but it also determines the volume of liquid 
that must flow across the heterogeneity and hence the flux of liquid through the rivulets that form 
as a result of depinning. Thus when L is large three-dimensional instabilities are always present 
and the depinning of a ridge is preceded by an instability of the 2d ridge that forms 3d drop-like 
states that remain pinned to the heterogeneity, with a further increase in the driving required before 
depinning takes place. 

Our results show that the bifurcation from 2d ridges to 3d drop-like states is subcritical, at least 
for the parameter values we use. As a result the 3d states are initially unstable, although with 
increasing amplitude the 3d states turn around towards larger values of fx and may acquire stability 
once they resemble steady rivulets. However, once the 2d ridges near the 2d sniper are unstable 
to 3d states, the 2d depinning mechanism no longer leads to stable stick-slip motion of the ridges, 
and instead we find a stable 3d version, which we call stick-slip drops. Our numerics suggest that 
in both the 2d and 3d cases the SSR and SSD states terminate in global bifurcations involving 
equilibria of saddle-type. However, the theory of such bifurcations shows that the periodic orbits 
involved in the global bifurcations are unstable and indeed our computations in the 2d case provide 
some evidence of a saddle-node bifurcation in the neighborhood of the putative global bifurcation 
at which the stable SSR states lose stability. Unfortunately we are unable to follow unstable time- 
periodic states such as the SSR and SSD states to confirm the presence of such global bifurcations. 

We were able, however, to interpret the transitions we observed using a two-dimensional model 
system analyzed by Krauskopf and Oldeman [901]. This system models the dynamics arising from 
the interaction of a sniper and a pitchfork, and its analysis reveals the great wealth of behavior 
available to systems of this type. We expect that transitions present for other values of our param- 
eters can likewise be interpreted in terms of this model system. Model systems such as Eqs. ©-([7]) 
suffer from a drawback, however, in that the dynamics is necessarily two-dimensional and hence 
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no chaos is present in the model. A three-dimensional extension of the model Ii92ll does exhibit 
chaotic dynamics associated with the plethora of global bifurcations revealed in [90] and offers a 
glimpse into the potential behavior associated with the (formally infinite-dimensional) depinning 
problem. However, our simulations have provided no conclusive evidence for chaotic dynamics. 
In some regimes we have also found a Hopf bifurcation in the 2d depinning problem, usually very 
close to both the saddle-node of the ridge states and to the 3d instability, suggesting that the full 
dynamics of the system can be captured by examining the interaction of the Takens-Bogdanov 



bifurcation with the pitchfork leading to 3d states Il93h . with reinjection as in [90]. 

In the parameter regime where the long-wave approximation applies the problem studied here, 
viz. film flow and drop motion on a heterogeneous substrate with a well-defined spatial period of 
the wettability, is closely related to the corresponding flow on the outside or inside of a rotating 
horizontal cylinder Ii52ll . In particular, the depinning dynamics of a two-dimensional drop of 
partially wetting liquid on the surface of a rotating cylinder has a close counterpart in the depinning 
dynamics via a sniper bifurcation described here for two-dimensional drops on heterogeneous 



substrates [|39l 



521] . Further exploration of the analogy between these two systems may therefore be 



fruitful. In particular, it may be possible to relate the three-dimensional structures and transitions 
described in the present paper to flows on rotating cylinders. For instance, the transition from 
stick-slip drops to spatially modulated rivulets via a Hopf bifurcation (Fig. [25]) would correspond 
to a transition from drops rotating with the cylinder (a state that has apparently not been described 
in the literature on film flow on the outside of a cylinder) to azimuthal rings [94]. 

Experiments on pattern formation in flows of thick films on the inner surface of a rotating 
cylinder report a variety of different structures and transitions that resemble some of our results 
|95,|9^]. These include, for instance, stationary straight and wavy fronts (similar to our steady 



spanwise invariant ridge and spanwise modulated ridge states). The latter are called "shark teeth" 
in J^i . At higher rotation velocities the stationary wavy fronts decay into drops that are "dragged 



up the receding wall" j95\\ . a transition corresponding to what we describe here as a 3d depinning 
transition. Stationary localized bumps, i.e., pinned drops, may also be present j95\ \. 
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